syms x theta p_H1 p_H0 p_L1 p_L0 theta_min theta_max
g_theta = 1 / (theta_max - theta_min);
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
disp('D_H1 = '); disp(D_H1)
disp('D_H0 = '); disp(D_H0)
disp('D_L1 = '); disp(D_L1)
disp('D_L0 = '); disp(D_L0)
% D_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% D_H0 = int(int(-exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% D_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% D_L0 = int(int(-exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
disp('Pi_H1 = '); disp(Pi_H1)
disp('Pi_H0 = '); disp(Pi_H0)
disp('Pi_L1 = '); disp(Pi_L1)
disp('Pi_L0 = '); disp(Pi_L0)
% Pi_H1 = -(c_H1 - p_H1)*int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% Pi_H0 = -int(int(-exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)*(c_H0 - p_H0);
% Pi_L1 = -(c_L1 - p_L1)*int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1);
% Pi_L0 = -int(int(-exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)*(c_L0 - p_L0);
% int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - int(int((lambda*exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1)*(c_H1 - p_H1) == 0
% int(int(-exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - (c_H0 - p_H0)*int(int((lambda*exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(- tau*x^2 + R - p_H0 + s_H*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1) == 0
% int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - int(int((lambda*exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1) == 0
% int(int(-exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))/
% ((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))
% + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta))
% + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))), theta, theta_min, theta_max), x, 0, 1)
% - (c_L0 - p_L0)*int(int((lambda*exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta))))
% - (lambda*exp(2*lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))/((theta_min - theta_max)
% *(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))
% + exp(lambda*(- tau*x^2 + R - p_H0 + s_H*theta)) + exp(lambda*(- tau*x^2 + R - p_L0 + s_L*theta)))^2)
% , theta, theta_min, theta_max), x, 0, 1) == 0
% to find the profit-maximizing prices p_H1*, p_H0*, p_L1*, and p_L0*. => define parameteres substitute into the equations and solve
syms p_H1 p_H0 p_L1 p_L0 theta x
% % Define parameter values
denominator = (theta_max - theta_min) * ( ...
exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + ...
exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + ...
exp(lambda*(-tau*x^2 + R - p_H0 + s_H*theta)) + ...
exp(lambda*(-tau*x^2 + R - p_L0 + s_L*theta)) );
foc_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_H1 - p_H1) == 0;
foc_H0 = int(int(-exp(lambda*(-tau*x^2 + R - p_H0 + s_H*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H0 - p_H0) * int(int((lambda * exp(lambda*(-tau*x^2 + R - p_H0 + s_H*theta))) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R - p_H0 + s_H*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
foc_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_L1 - p_L1) == 0;
foc_L0 = int(int(-exp(lambda*(-tau*x^2 + R - p_L0 + s_L*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L0 - p_L0) * int(int((lambda * exp(lambda*(-tau*x^2 + R - p_L0 + s_L*theta))) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R - p_L0 + s_L*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
solution = vpasolve([foc_H1, foc_H0, foc_L1, foc_L0], [p_H1, p_H0, p_L1, p_L0]);
disp(['p_H1* = ', char(solution.p_H1)]);
p_H1* = 0.97818251806516813907110255397178
disp(['p_H0* = ', char(solution.p_H0)]);
p_H0* = 0.97818251806516813907110255397178
disp(['p_L1* = ', char(solution.p_L1)]);
p_L1* = 0.67299421682353783913406045511087
disp(['p_L0* = ', char(solution.p_L0)]);
p_L0* = 0.67299421682353783913406045511087
denominator = @(x, theta) (theta_max - theta_min) .* ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) );
D_H1 = integral2(@(theta, x) exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
fprintf('D_H1 = %.4f\n', D_H1);
fprintf('D_H0 = %.4f\n', D_H0);
fprintf('D_L1 = %.4f\n', D_L1);
fprintf('D_L0 = %.4f\n', D_L0);
fprintf('Pi_H1 = %.4f\n', Pi_H1);
fprintf('Pi_H0 = %.4f\n', Pi_H0);
fprintf('Pi_L1 = %.4f\n', Pi_L1);
fprintf('Pi_L0 = %.4f\n', Pi_L0);
p_H1 = 0.97; p_H0 = 0.97; p_L1 = 0.67; p_L0 = 0.67;
theta_samples = 0.2 + (2.2 - 0.2) * rand(N,1);
U_H1 = R + s_H * theta_samples - p_H1 - tau * (1 - x_samples).^2;
U_H0 = R + s_H * theta_samples - p_H0 - tau * x_samples.^2;
U_L1 = R + s_L * theta_samples - p_L1 - tau * (1 - x_samples).^2;
U_L0 = R + s_L * theta_samples - p_L0 - tau * x_samples.^2;
denominator = exp(lambda * U_H1) + exp(lambda * U_H0) + exp(lambda * U_L1) + exp(lambda * U_L0);
P_H1 = exp(lambda * U_H1) ./ denominator;
P_H0 = exp(lambda * U_H0) ./ denominator;
P_L1 = exp(lambda * U_L1) ./ denominator;
P_L0 = exp(lambda * U_L0) ./ denominator;
scatter(x_samples, theta_samples, 10, P_H1, 'filled');
title('Probability for H1');
scatter(x_samples, theta_samples, 10, P_H0, 'filled');
title('Probability for H0');
scatter(x_samples, theta_samples, 10, P_L1, 'filled');
title('Probability for L1');
scatter(x_samples, theta_samples, 10, P_L0, 'filled');
title('Probability for L0');
syms x theta p_H1 p_H0 p_L1 p_L0 theta_min theta_max
g_theta = 1 / (theta_max - theta_min);
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
disp('D_H1 = '); disp(D_H1)
disp('D_H0 = '); disp(D_H0)
disp('D_L1 = '); disp(D_L1)
disp('D_L0 = '); disp(D_L0)
% D_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% D_H0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% D_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% D_L0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% tax just applied in demand!
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
disp('Pi_H1 = '); disp(Pi_H1)
disp('Pi_H0 = '); disp(Pi_H0)
disp('Pi_L1 = '); disp(Pi_L1)
disp('Pi_L0 = '); disp(Pi_L0)
% Pi_H1 = -int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_H1 - p_H1);
% Pi_H0 = -(c_H0 - p_H0)*int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% Pi_L1 = -int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1);
% Pi_L0 = -(c_L0 - p_L0)*int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
dPi_H1_dpH1 = diff(Pi_H1, p_H1) == 0;
dPi_H0_dpH0 = diff(Pi_H0, p_H0) == 0;
dPi_L1_dpL1 = diff(Pi_L1, p_L1) == 0;
dPi_L0_dpL0 = diff(Pi_L0, p_L0) == 0;
disp('dPi_H1_dpH1:'); disp(dPi_H1_dpH1)
disp('dPi_H0_dpH0:'); disp(dPi_H0_dpH0)
disp('dPi_L1_dpL1:'); disp(dPi_L1_dpL1)
disp('dPi_L0_dpL0:'); disp(dPi_L0_dpL0)
% dPi_H1_dpH1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_H1 - p_H1)*int(int((lambda*exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_H0_dpH0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - int(int((lambda*exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1)*(c_H0 - p_H0) == 0;
% dPi_L1_dpL1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_L1 - p_L1)*int(int((lambda*exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_L0_dpL0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - int(int((lambda*exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1)*(c_L0 - p_L0) == 0;
syms p_H1 p_H0 p_L1 p_L0 theta x
denominator = (theta_max - theta_min) * ( ...
exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) + ...
exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + ...
exp(lambda*(-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + ...
exp(lambda*(-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) );
foc_H1 = int(int(-exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H1 - p_H1) * int(int((lambda * exp(lambda*(R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_H1 - tau*(x - 1)^2 + s_H*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
foc_H0 = int(int(-exp(lambda*(-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R + s_H*theta - p_H0*(t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_H0 - p_H0) == 0;
foc_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L1 - p_L1) * int(int((lambda * exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_L1 - tau*(x - 1)^2 + s_L*theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
foc_L0 = int(int(-exp(lambda*(-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- int(int((lambda * exp(lambda*(-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau*x^2 + R + s_L*theta - p_L0*(t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) * (c_L0 - p_L0) == 0;
solution = vpasolve([foc_H1, foc_H0, foc_L1, foc_L0], [p_H1, p_H0, p_L1, p_L0]);
disp(['p_H1** = ', char(solution.p_H1)]);
p_H1** = 0.98245738528891675184015215351829
disp(['p_H0** = ', char(solution.p_H0)]);
p_H0** = 0.8451818220778982201818734232353
disp(['p_L1** = ', char(solution.p_L1)]);
p_L1** = 0.67507753269566794239079645081273
disp(['p_L0** = ', char(solution.p_L0)]);
p_L0** = 0.5602969265225777050338138212517
denominator = @(x, theta) (theta_max - theta_min) .* ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0*(1+t))) + ...
exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0*(1+t))) );
D_H1 = integral2(@(theta, x) exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0*(1+t))) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0*(1+t))) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
fprintf('D_H1 = %.4f\n', D_H1);
fprintf('D_H0 = %.4f\n', D_H0);
fprintf('D_L1 = %.4f\n', D_L1);
fprintf('D_L0 = %.4f\n', D_L0);
fprintf('Pi_H1 = %.4f\n', Pi_H1);
fprintf('Pi_H0 = %.4f\n', Pi_H0);
fprintf('Pi_L1 = %.4f\n', Pi_L1);
fprintf('Pi_L0 = %.4f\n', Pi_L0);
p_H0_tax = p_H0 * (1 + t);
p_L0_tax = p_L0 * (1 + t);
fprintf('\nConsumer Prices:\n');
fprintf('p_H0_VAT = %.4f\n', p_H0_tax);
fprintf('p_L0_VAT = %.4f\n', p_L0_tax);
syms x theta p_H1 p_H0 p_L1 p_L0 theta_min theta_max
g_theta = 1 / (theta_max - theta_min);
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
disp('D_H1 = '); disp(D_H1)
disp('D_H0 = '); disp(D_H0)
disp('D_L1 = '); disp(D_L1)
disp('D_L0 = '); disp(D_L0)
% D_H1 = int(int(-exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1
% D_H0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)
% D_L1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)
% D_L0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
disp('Pi_H1 = '); disp(Pi_H1)
disp('Pi_H0 = '); disp(Pi_H0)
disp('Pi_L1 = '); disp(Pi_L1)
disp('Pi_L0 = '); disp(Pi_L0)
% Pi_H1 = -(c_H1 - p_H1)*int(int(-exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1);
% Pi_H0 = -int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_H0 - p_H0);
% Pi_L1 = -int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1);
% Pi_L0 = -int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1)*(c_L0 - p_L0);
dPi_H1_dpH1 = diff(Pi_H1, p_H1) == 0;
dPi_H0_dpH0 = diff(Pi_H0, p_H0) == 0;
dPi_L1_dpL1 = diff(Pi_L1, p_L1) == 0;
dPi_L0_dpL0 = diff(Pi_L0, p_L0) == 0;
disp('dPi_H1_dpH1 = '); disp(dPi_H1_dpH1)
disp('dPi_H0_dpH0 = '); disp(dPi_H0_dpH0)
disp('dPi_L1_dpL1 = '); disp(dPi_L1_dpL1)
disp('dPi_L0_dpL0 = '); disp(dPi_L0_dpL0)
% dPi_H1_dpH1 = int(int(-exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_H1 - p_H1)*int(int((lambda*exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_H0_dpH0 = int(int(-exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_H0 - p_H0)*int(int((lambda*exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
% dPi_L1_dpL1 = int(int(-exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - int(int((lambda*exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1)*(c_L1 - p_L1) == 0;
% dPi_L0_dpL0 = int(int(-exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))), theta, theta_min, theta_max), x, 0, 1) - (c_L0 - p_L0)*int(int((lambda*exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))) - (lambda*exp(2*lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1)))*(t + 1))/((theta_min - theta_max)*(exp(lambda*(R - p_L1 - tau*(x - 1)^2 + s_L*theta)) + exp(lambda*(- tau*x^2 + R + s_H*theta - p_H0*(t + 1))) + exp(lambda*(- tau*x^2 + R + s_L*theta - p_L0*(t + 1))) + exp(lambda*(R - tau*(x - 1)^2 + s_H*theta - p_H1*(t + 1))))^2), theta, theta_min, theta_max), x, 0, 1) == 0;
syms p_H1 p_H0 p_L1 p_L0 theta x
denominator = (theta_max - theta_min) * ( ...
exp(lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta)) + ...
exp(lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) + ...
exp(lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) + ...
exp(lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) );
foc_H1 = int(int(-exp(lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H1 - p_H1) * int(int((lambda * exp(lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (R - tau * (x - 1)^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
foc_H0 = int(int(-exp(lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_H0 - p_H0) * int(int((lambda * exp(lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau * x^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
foc_L1 = int(int(-exp(lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta)) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L1 - p_L1) * int(int((lambda * exp(lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta))) / denominator ...
- (lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1)^2 + s_L * theta))) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
foc_L0 = int(int(-exp(lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) / denominator, ...
theta, theta_min, theta_max), x, 0, 1) ...
- (c_L0 - p_L0) * int(int((lambda * exp(lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) / denominator ...
- (lambda * exp(2 * lambda * (-tau * x^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) / denominator^2, ...
theta, theta_min, theta_max), x, 0, 1) == 0;
solution = vpasolve([foc_H1, foc_H0, foc_L1, foc_L0], [p_H1, p_H0, p_L1, p_L0]);
disp(['p_H1*** = ', char(solution.p_H1)]);
p_H1*** = 0.84730365281426619546342786138934
disp(['p_H0*** = ', char(solution.p_H0)]);
p_H0*** = 0.84809017803927929379959703727505
disp(['p_L1*** = ', char(solution.p_L1)]);
p_L1*** = 0.67999498272950621551040317326515
disp(['p_L0*** = ', char(solution.p_L0)]);
p_L0*** = 0.56166362780346060630817744594582
p_H1_tax = p_H1 * (1 + t);
p_H0_tax = p_H0 * (1 + t);
p_L0_tax = p_L0 * (1 + t);
denominator = @(x, theta) (theta_max - theta_min) .* ( ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0_tax)) + ...
exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0_tax)) + ...
exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1_tax)));
D_H1 = integral2(@(theta, x) exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1_tax)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0_tax)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0_tax)) ./ denominator(x, theta), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
Pi_H1 = (p_H1 - c_H1) * D_H1;
Pi_H0 = (p_H0 - c_H0) * D_H0;
Pi_L1 = (p_L1 - c_L1) * D_L1;
Pi_L0 = (p_L0 - c_L0) * D_L0;
fprintf('D_H1 = %.4f\n', D_H1);
fprintf('D_H0 = %.4f\n', D_H0);
fprintf('D_L1 = %.4f\n', D_L1);
fprintf('D_L0 = %.4f\n', D_L0);
fprintf('Pi_H1 = %.4f\n', Pi_H1);
fprintf('Pi_H0 = %.4f\n', Pi_H0);
fprintf('Pi_L1 = %.4f\n', Pi_L1);
fprintf('Pi_L0 = %.4f\n', Pi_L0);
fprintf('\nConsumer Prices:\n');
fprintf('p_H1_VAT = %.4f\n', p_H1_tax);
fprintf('p_H0_VAT = %.4f\n', p_H0_tax);
fprintf('p_L0_VAT = %.4f\n', p_L0_tax);
theta_90 = theta_min + 0.9 * (theta_max - theta_min);
theta_10 = theta_min + 0.1 * (theta_max - theta_min);
denominator = @(theta,x) (exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) + exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) );
D_H1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) ./ denominator(theta,x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0)) ./ denominator(theta,x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
fprintf('D_H1 (Top 10%% WTP) = %.4f\n', D_H1_top10);
D_H1 (Top 10% WTP) = 0.2040
fprintf('D_H0 (Top 10%% WTP) = %.4f\n', D_H0_top10);
D_H0 (Top 10% WTP) = 0.2040
fprintf('D_L1 (Top 10%% WTP) = %.4f\n', D_L1_top10);
D_L1 (Top 10% WTP) = 0.0460
fprintf('D_L0 (Top 10%% WTP) = %.4f\n', D_L0_top10);
D_L0 (Top 10% WTP) = 0.0460
fprintf('D_H1 (Bottom 10%% WTP) = %.4f\n', D_H1_low10);
D_H1 (Bottom 10% WTP) = 0.1058
fprintf('D_H0 (Bottom 10%% WTP) = %.4f\n', D_H0_low10);
D_H0 (Bottom 10% WTP) = 0.1058
fprintf('D_L1 (Bottom 10%% WTP) = %.4f\n', D_L1_low10);
D_L1 (Bottom 10% WTP) = 0.1442
fprintf('D_L0 (Bottom 10%% WTP) = %.4f\n', D_L0_low10);
D_L0 (Bottom 10% WTP) = 0.1442
theta_90 = theta_min + 0.9 * (theta_max - theta_min);
theta_10 = theta_min + 0.1 * (theta_max - theta_min);
denominator = @(theta, x) ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) + ...
exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) );
D_H1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
fprintf('D_H1 (Top 10%% WTP) = %.4f\n', D_H1_top10);
D_H1 (Top 10% WTP) = 0.2138
fprintf('D_H0 (Top 10%% WTP) = %.4f\n', D_H0_top10);
D_H0 (Top 10% WTP) = 0.1902
fprintf('D_L1 (Top 10%% WTP) = %.4f\n', D_L1_top10);
D_L1 (Top 10% WTP) = 0.0485
fprintf('D_L0 (Top 10%% WTP) = %.4f\n', D_L0_top10);
D_L0 (Top 10% WTP) = 0.0475
fprintf('D_H1 (Bottom 10%% WTP) = %.4f\n', D_H1_low10);
D_H1 (Bottom 10% WTP) = 0.1091
fprintf('D_H0 (Bottom 10%% WTP) = %.4f\n', D_H0_low10);
D_H0 (Bottom 10% WTP) = 0.0962
fprintf('D_L1 (Bottom 10%% WTP) = %.4f\n', D_L1_low10);
D_L1 (Bottom 10% WTP) = 0.1495
fprintf('D_L0 (Bottom 10%% WTP) = %.4f\n', D_L0_low10);
D_L0 (Bottom 10% WTP) = 0.1452
theta_90 = theta_min + 0.9 * (theta_max - theta_min);
theta_10 = theta_min + 0.1 * (theta_max - theta_min);
denominator = @(theta, x) (exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (1 + t))) + exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) + exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) );
D_H1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 * (1 + t) - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0_top10 = (1 / (theta_max - theta_90)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_90, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_H1 * (1 + t) - tau * (x - 1).^2 + s_H * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0_low10 = (1 / (theta_10 - theta_min)) * integral2(@(theta, x) ...
(1 / (theta_max - theta_min)) * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (1 + t))) ./ denominator(theta, x), ...
theta_min, theta_10, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
fprintf('D_H1 (Top 10%% WTP) = %.4f\n', D_H1_top10);
D_H1 (Top 10% WTP) = 0.1989
fprintf('D_H0 (Top 10%% WTP) = %.4f\n', D_H0_top10);
D_H0 (Top 10% WTP) = 0.1991
fprintf('D_L1 (Top 10%% WTP) = %.4f\n', D_L1_top10);
D_L1 (Top 10% WTP) = 0.0521
fprintf('D_L0 (Top 10%% WTP) = %.4f\n', D_L0_top10);
D_L0 (Top 10% WTP) = 0.0499
fprintf('D_H1 (Bottom 10%% WTP) = %.4f\n', D_H1_low10);
D_H1 (Bottom 10% WTP) = 0.0979
fprintf('D_H0 (Bottom 10%% WTP) = %.4f\n', D_H0_low10);
D_H0 (Bottom 10% WTP) = 0.0983
fprintf('D_L1 (Bottom 10%% WTP) = %.4f\n', D_L1_low10);
D_L1 (Bottom 10% WTP) = 0.1547
fprintf('D_L0 (Bottom 10%% WTP) = %.4f\n', D_L0_low10);
D_L0 (Bottom 10% WTP) = 0.1490
syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda theta_min theta_max
U_H1 = R + theta*s_H - p_H1 - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0 - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0 - tau*x^2;
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
g_theta = 1 / (theta_max - theta_min);
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
disp('D_H1 = '); disp(D_H1)
disp('D_H0 = '); disp(D_H0)
disp('D_L1 = '); disp(D_L1)
disp('D_L0 = '); disp(D_L0)
dD_H1_dp_H1 = diff(D_H1, p_H1);
dD_H0_dp_H0 = diff(D_H0, p_H0);
dD_L1_dp_L1 = diff(D_L1, p_L1);
dD_L0_dp_L0 = diff(D_L0, p_L0);
e_H1 = simplify(dD_H1_dp_H1 * (p_H1/D_H1));
e_H0 = simplify(dD_H0_dp_H0 * (p_H0/D_H0));
e_L1 = simplify(dD_L1_dp_L1 * (p_L1/D_L1));
e_L0 = simplify(dD_L0_dp_L0 * (p_L0/D_L0));
disp('e_H1 ='); disp(e_H1)
disp('e_H0 ='); disp(e_H0)
disp('e_L1 ='); disp(e_L1)
disp('e_L0 ='); disp(e_L0)
denominator = @(x, theta) ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) + ...
exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) );
D_H1 = integral2(@(theta, x) - exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0 = integral2(@(theta, x) exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_H1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_H0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R - p_H0 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R - p_H0 + s_H * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_L1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_L0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R - p_L0 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R - p_L0 + s_L * theta))) ./ ((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
e_H1 = dD_H1_dP * (p_H1 / D_H1);
e_H0 = dD_H0_dP * (p_H0 / D_H0);
e_L1 = dD_L1_dP * (p_L1 / D_L1);
e_L0 = dD_L0_dP * (p_L0 / D_L0);
disp('Computed Demand Values:')
fprintf('D_H1 = %.6f\n', D_H1);
fprintf('D_H0 = %.6f\n', D_H0);
fprintf('D_L1 = %.6f\n', D_L1);
fprintf('D_L0 = %.6f\n', D_L0);
disp('Computed Demand Derivatives:')
Computed Demand Derivatives:
fprintf('dD_H1/dP = %.6f\n', dD_H1_dP);
fprintf('dD_H0/dP = %.6f\n', dD_H0_dP);
fprintf('dD_L1/dP = %.6f\n', dD_L1_dP);
fprintf('dD_L0/dP = %.6f\n', dD_L0_dP);
disp('Computed Elasticities:')
fprintf('e_H1 = %.4f\n', e_H1);
fprintf('e_H0 = %.4f\n', e_H0);
fprintf('e_L1 = %.4f\n', e_L1);
fprintf('e_L0 = %.4f\n', e_L0);
syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda theta_min theta_max t
U_H1 = R + theta*s_H - p_H1 - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0*(1+t) - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0*(1+t) - tau*x^2;
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
g_theta = 1 / (theta_max - theta_min);
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
dD_H1_dp_H1 = diff(D_H1, p_H1);
dD_H0_dp_H0 = diff(D_H0, p_H0);
dD_L1_dp_L1 = diff(D_L1, p_L1);
dD_L0_dp_L0 = diff(D_L0, p_L0);
e_H1 = simplify(dD_H1_dp_H1 * (p_H1/D_H1));
e_H0 = simplify(dD_H0_dp_H0 * (p_H0/D_H0));
e_L1 = simplify(dD_L1_dp_L1 * (p_L1/D_L1));
e_L0 = simplify(dD_L0_dp_L0 * (p_L0/D_L0));
disp('e_H1 ='); disp(e_H1)
disp('e_H0 ='); disp(e_H0)
disp('e_L1 ='); disp(e_L1)
disp('e_L0 ='); disp(e_L0)
denominator = @(x, theta) ( ...
exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (R - tau * x.^2 + s_H * theta - (1 + t) * p_H0)) + ...
exp(lambda * (R - tau * x.^2 + s_L * theta - (1 + t) * p_L0)) );
D_H1 = integral2(@(theta, x) exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0 = integral2(@(theta, x) exp(lambda * (R - tau * x.^2 + s_H * theta - (1 + t) * p_H0)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1 = integral2(@(theta, x) exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0 = integral2(@(theta, x) exp(lambda * (R - tau * x.^2 + s_L * theta - (1 + t) * p_L0)) ./ ...
((theta_max - theta_min) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_H1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_H1 - tau * (x - 1).^2 + s_H * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_H0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_L1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_L0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
e_H1 = dD_H1_dP * (p_H1 / D_H1);
e_H0 = dD_H0_dP * (p_H0 / D_H0);
e_L1 = dD_L1_dP * (p_L1 / D_L1);
e_L0 = dD_L0_dP * (p_L0 / D_L0);
disp('Computed Demand Values with Tax:')
Computed Demand Values with Tax:
fprintf('D_H1 = %.6f\n', D_H1);
fprintf('D_H0 = %.6f\n', D_H0);
fprintf('D_L1 = %.6f\n', D_L1);
fprintf('D_L0 = %.6f\n', D_L0);
disp('Computed Demand Derivatives with Tax:')
Computed Demand Derivatives with Tax:
fprintf('dD_H1/dP = %.6f\n', dD_H1_dP);
fprintf('dD_H0/dP = %.6f\n', dD_H0_dP);
fprintf('dD_L1/dP = %.6f\n', dD_L1_dP);
fprintf('dD_L0/dP = %.6f\n', dD_L0_dP);
disp('Computed Elasticities with Tax:')
Computed Elasticities with Tax:
fprintf('e_H1 = %.4f\n', e_H1);
fprintf('e_H0 = %.4f\n', e_H0);
fprintf('e_L1 = %.4f\n', e_L1);
fprintf('e_L0 = %.4f\n', e_L0);
syms x theta p_H1 p_H0 p_L1 p_L0 tau R s_H s_L lambda theta_min theta_max t
U_H1 = R + theta*s_H - p_H1*(1+t) - tau*(1 - x)^2;
U_H0 = R + theta*s_H - p_H0*(1+t) - tau*x^2;
U_L1 = R + theta*s_L - p_L1 - tau*(1 - x)^2;
U_L0 = R + theta*s_L - p_L0*(1+t) - tau*x^2;
denominator = exp(lambda*U_H1) + exp(lambda*U_H0) + exp(lambda*U_L1) + exp(lambda*U_L0);
Pr_H1 = exp(lambda*U_H1) / denominator;
Pr_H0 = exp(lambda*U_H0) / denominator;
Pr_L1 = exp(lambda*U_L1) / denominator;
Pr_L0 = exp(lambda*U_L0) / denominator;
g_theta = 1 / (theta_max - theta_min);
D_H1 = int(int(Pr_H1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_H0 = int(int(Pr_H0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L1 = int(int(Pr_L1 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
D_L0 = int(int(Pr_L0 * f_x * g_theta, theta, theta_min, theta_max), x, 0, 1);
disp('D_H1 = '); disp(D_H1)
disp('D_H0 = '); disp(D_H0)
disp('D_L1 = '); disp(D_L1)
disp('D_L0 = '); disp(D_L0)
dD_H1_dp_H1 = diff(D_H1, p_H1);
dD_H0_dp_H0 = diff(D_H0, p_H0);
dD_L1_dp_L1 = diff(D_L1, p_L1);
dD_L0_dp_L0 = diff(D_L0, p_L0);
e_H1 = simplify(dD_H1_dp_H1 * (p_H1/D_H1));
e_H0 = simplify(dD_H0_dp_H0 * (p_H0/D_H0));
e_L1 = simplify(dD_L1_dp_L1 * (p_L1/D_L1));
e_L0 = simplify(dD_L0_dp_L0 * (p_L0/D_L0));
disp('e_H1 ='); disp(e_H1)
disp('e_H0 ='); disp(e_H0)
disp('e_L1 ='); disp(e_L1)
disp('e_L0 ='); disp(e_L0)
denominator = @(x, theta) ( ...
exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - (1 + t) * p_H1)) + ...
exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) + ...
exp(lambda * (R - tau * x.^2 + s_H * theta - (1 + t) * p_H0)) + ...
exp(lambda * (R - tau * x.^2 + s_L * theta - (1 + t) * p_L0)) );
D_H1 = integral2(@(theta, x) ...
-exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (t + 1))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_H0 = integral2(@(theta, x) ...
-exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L1 = integral2(@(theta, x) ...
-exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
D_L0 = integral2(@(theta, x) ...
-exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)), ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_H1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - tau * (x - 1).^2 + s_H * theta - p_H1 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_H0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_H * theta - p_H0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_L1_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (R - p_L1 - tau * (x - 1).^2 + s_L * theta))) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
dD_L0_dP = integral2(@(theta, x) ...
(lambda * exp(lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)) - ...
(lambda * exp(2 * lambda * (-tau * x.^2 + R + s_L * theta - p_L0 * (t + 1))) * (t + 1)) ./ ...
((theta_min - theta_max) .* denominator(x, theta)).^2, ...
theta_min, theta_max, 0, 1, 'RelTol', 1e-6, 'AbsTol', 1e-6);
e_H1 = dD_H1_dP * (p_H1 / D_H1);
e_H0 = dD_H0_dP * (p_H0 / D_H0);
e_L1 = dD_L1_dP * (p_L1 / D_L1);
e_L0 = dD_L0_dP * (p_L0 / D_L0);
disp('Computed Demand Values with Tax:')
Computed Demand Values with Tax:
fprintf('D_H1 = %.6f\n', D_H1);
fprintf('D_H0 = %.6f\n', D_H0);
fprintf('D_L1 = %.6f\n', D_L1);
fprintf('D_L0 = %.6f\n', D_L0);
disp('Computed Demand Derivatives with Tax:')
Computed Demand Derivatives with Tax:
fprintf('dD_H1/dP = %.6f\n', dD_H1_dP);
fprintf('dD_H0/dP = %.6f\n', dD_H0_dP);
fprintf('dD_L1/dP = %.6f\n', dD_L1_dP);
fprintf('dD_L0/dP = %.6f\n', dD_L0_dP);
disp('Computed Elasticities with Tax:')
Computed Elasticities with Tax:
fprintf('e_H1 = %.4f\n', e_H1);
fprintf('e_H0 = %.4f\n', e_H0);
fprintf('e_L1 = %.4f\n', e_L1);
fprintf('e_L0 = %.4f\n', e_L0);